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ADVECTION-DIFFUSION EQUATIONS 
WITH DENSITY CONSTRAINTS 

ALPAr RICHARD MESZAROS AND FILIPPO SANTAMBROGIO 


Abstract. In the spirit of the macroscopic crowd motion models with hard congestion (i.e. a 
strong density constraint p < 1) introduced by Maury et al. some years ago, we analyze a variant 
of the same models where diffusion of the agents is also taken into account. From the modeling 
point of view, this means that individuals try to follow a given spontaneous velocity, but are subject 
to a Brownian diffusion, and have to adapt to a density constraint which introduces a pressure term 
affecting the movement. From the PDE point of view, this corresponds to a modified Fokker-Planck 
equation, with an additional gradient of a pressure (only living in the saturated zone {p = 1}) in 
the drift. The paper proves existence and some estimates, based on optimal transport techniques. 


1. Introduction 

In the past few years modeling crowd behavior has become a very active field of applied mathe¬ 
matics. Beyond their importance in real life applications, these modeling problems serve as basic 
ideas to understand many other phenomena coming for example from biology (cell migration, tumor 
growth, pattern formations in animal populations, etc.), particle physics and economics. A first 
non-exhaustive list of references for these problems is [9, 11, 12, 13, 17, 19, 20, 22, 23, 32], A very 
natural question in all these models is the problem of congestion phenomenon: in many practical 
situations, very high quantities of individuals could try to occupy the same spot, which could be 
impossible, or lead to strong negative effects on the motion, because of natural limitations on the 
crowd density. 

These phenomena have been studied by using different models, which could be either “micro¬ 
scopic” (based on ODEs on the motion of a high number of agents) or “macroscopic” (describing 
the agents via their density and velocity, typically with Eulerian formalism). Let us concentrate on 
the macroscopic models, where the density p plays a crucial role. These very same models can be 
characterized either by “soft congestion” effects (i.e. the higher the density the slower the motion), 
or by “hard congestion” (i.e. an abrupt threshold effect: if the density touches a certain maximal 
value, the motion is strongly affected, while nothing happens for smaller values of the density). 
See [31] for comparison between the different classes of models. This last class of models, due to 
the discontinuity in the congestion effects, presents new mathematical difficulties, which cannot be 
analyzed with the usual techniques from conservation laws (or, more generally, evolution PDEs) 
used for soft congestion. 

A very powerful tool to attack macroscopic hard-congestion problems is the theory of optimal 
transportation (see [41, 40]), as we can see in [30, 31, 37, 38]. In this framework, the density of the 
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agents solves a continuity equation (with velocity field taking into account the congestion effects), 
and can be seen as a curve in the Wasserstein space. 

Our aim in this paper is to endow the macroscopic hard congestion models of [30, 31, 37, 38] 
with diffusion effects. In other words, we will study an evolution equation where particles 

• have a spontaneous velocity field ut(x) which depends on time and on their position, and 
is the velocity they would follow in the absence of the other particles, 

• must adapt their velocity to the existence of an incompressibility constraint which prevents 
the density to go beyond a given threshold, 

• are subject to some diffusion effect. 

This can be considered as a model for a crowd where a part of the motion of each agent is driven 
by a Brownian motion. Implementing this new element into the existing models could give a better 
approximation of reality: as usual when one adds a stochastic component, this can be a (very) 
rough approximation of unpredictable effects which are not already handled by the model, and this 
could work well when dealing with large populations. 

Anyway, we do not want to discuss here the validity of this hard-congestion model and we are 
mainly concerned with its mathematical analysis. In particular, we will consider existence and 
regularity estimates, while we do not treat the uniqueness issue. Uniqueness is considered in a 
recent work of the first author in collaboration with S. Di Marino, see [16], and one can observe 
that the insertion of diffusion dramatically simplifies the picture as far as uniqueness is concerned. 

We also underline that one of the goals of the current paper (and of [16]) is to better “prepare” 
these hard congestion crowd motion models for a possible analysis in the framework of Mean Field 
Games (see [26, 27, 28], and also [39]). These MFG models usually involve a stochastic term, also 
implying regularizing effects, that are useful in the mathematical analysis of the corresponding 
PDEs. 

1.1. The existing first order models in the light of [30, 31]. Some macroscopic models for 
crowd motion with density constraints and “hard congestion” effects were studied in [31] and [30]. 
We briefly present them as follows: 

• The density of the population in a bounded (convex) domain Q C M. d is described by 
a probability measure p £ V(Q). The initial density po £ V(£ 1) evolves in time, and pt 
denotes its value at each time t € [0, T]. 

• The spontaneous velocity field of the population is a given time-dependent field, denoted 
by ut■ It represents the velocity that each individual would like to follow in the absence of 
the others. Ignoring the density constraint, this would give rise to the continuity equation 
dtPt + V • (pm) = 0. We observe that in the original work [30] the vector field ut(x) was 
taken of the form -VD(x) (independent of time and of gradient form) but we try here 
to be more general (see [37] where the non-gradient case is studied under some stronger 
regularity assumptions). 

• The set of admissible densities will be denoted by K, := {p £ T’(fl) : p < 1}. In order to 
guarantee that K, is neither empty nor trivial, we suppose |U| > 1. 

• The set of admissible velocity fields with respect to the density p is characterized by the 
sign of the divergence of the velocity held on the saturated zone. We need to suppose also 
that all admissible velocity fields are such that no mass exists from the domain. So formally 
we set 

adrn(p) := ju : Q —» BA* : V • v > 0 on {p = 1} and v ■ n < 0 on <9f2 j . 
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We consider the projection operator P in L 2 {C d ): 


P, 


adm(p) Utl ^ al § mln i)eadm(p) 


\U — V 


2 dx. 


( 1 . 1 ) 


Note that we could have used the Hilbert space L 2 (p) instead of L 2 {C d ): this would be more 
natural in this kind of evolution equations, as L 2 (p ) is interpreted in a standard way as the 
tangent space to the Wasserstein space W2(17). Yet, these two projections turn out to be 
the same in this case, as the only relevant zone is {p = 1}. This is just formal, and would 
require more rigorous definitions (in particular of the divergence constraint in adm(p), see 
below). Anyway, to clarify, we choose to use the L 2 (£ d )-projection: in this way the vector 
fields are considered as defined Lebesgue-a.e. on the whole 17 (and not only on {p > 0}) 
and the dependence of the projected vector field on p only passes through the set adm(p). 
Finally we solve the following modified continuity equation for p 

9f Pt T V ' (Pl^ > adm(pt) [Vt]) 0, 


where the main point is that p is advected by a vector field, compatible with the constraints, 
which is the closest to the spontaneous one. 

The problem in solving Equation (1.1) is that the projected field has very low regularity: it is a 
priori only L 2 in x, and it does not depend smoothly on p either (since a density 1 and a density 
1 — e give very different projection operators). By the way, its divergence is not well-defined either. 
To handle this issue we need to redefine the set of admissible velocities by duality. Taking a test 
function p E H l (Q), p > 0 a.e., we obtain by the integration-by-parts equality 


v-Vpdx = — / (V-v)pdx + / pv ■ n dT~i d 1 (x). 

! Jfl J dfl 


For vector fields v which do not let mass go through the boundary <917 we have (in an a.e. sense) 
v ■ n = 0. This leads to the following definition 

adm(p) = E L 2 (17; R d ) : J v ■ Vp dx < 0, Vp € H l (p),p > 0,_p(l — p) = 0 a.e.| , 

(indeed, for smooth vector held with vanishing normal component on the boundary, this is equiv¬ 
alent to imposing V • v > 0 on the set {p = 1}). 

Now, if we set 

press(p) := {p € F7 1 (l7) : p > 0, p{l — p) = 0 a.e.} , 

we observe that, by definition, adm(p) and Vpress(p) are two convex cones which are dual to each 
other in L 2 (!7;R d ). Hence we always have a unique orthogonal decomposition 


( 1 . 2 ) 


u = v + Vp, v E adm(p); p E press(p), 


v ■ Vp dx = 0. 


In this decomposition (as it is the case every time we decompose on two dual convex cones), 
v = P a <im(p) ['“]■ These will be our mathematical definitions for adm(p) and for the projection onto 
this cone. 

Via this approach (introducing the new variable p and using its characterization from the previous 
line), for a given desired velocity held u : [0,T] x 17 —>• M d , the continuity equation (1.1) can be 
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rewritten as a system for the pair of variables (p,p) which is 

( d t pt + V • (p t (ut - Vp t )) = 0, in [0, T] x Cl, 

(1.3) < p > 0, p < 1, p{ 1 — p) = 0, in [0, T] x Cl, 

{ p t (u t - Vp t ) ■ n = 0, on [0, T] x dCl. 

This system is endowed with the initial condition p(0,x) = po(x) (po S 1C). As far as the spatial 
boundary dCl is concerned, we put no-flux boundary conditions to preserve the mass in Cl. 

Note that in the above system we withdrew the condition f (ut — V pt ) • Xpt = 0, as it is a 
consequence of the system (1.3) itself. Informally, this can be seen in the following way: for an 
arbitrary po € press (pt 0 ), we have that t i-A J n poPt is maximal at t = to (where it is equal to 
JqPo)- Differentiating this quantity w.r.t. t at t = to, using the equation (1.3), we get the desired 
orthogonality condition at t = to. For a rigorous proof of this fact (which holds for a.e. to), we 
refer to Proposition 4.7 in [15]. 


1.2. A diffusive counterpart. The goal of our work is to study a second order model of crowd 
movements with hard congestion effects where beside the transport factor a non-degenerate diffusion 
is present as well. The diffusion is the consequence of a randomness (a Brownian motion) in the 
movement of the crowd. 

With the ingredients that we introduced so far, we will modify the Fokker-Planck equation 
dtpt — A p t + V • ( ptUt) = 0 in order to take into account the density constraint pt < 1. Assuming 
enough regularity for the velocity field u, we observe that the Fokker Planck equation is derived from 
a motion given by the SODE dX t = ut(X t ) dt + \/2d B t (where B t is the standard d-dimensional 
Brownian motion), but is macroscopically represented by the advection of the density pt by the 
vector held — Vpt/pt + ut■ Projecting onto the set of admissible velocities raises a natural question: 
should we project only ut, and then apply the diffusion, or project the whole vector held, including 
—Vpt/pt? But this is not a real issue, since, at least formally, V pt/pt = 0 on the saturated 
Set { Pl, 1} and P&dm(pt)[ ^Pt/pt T U/J ^adm(pt)[ ^ Pt/pt\ T 0 T ■ 

Rigorously, this corresponds to the fact that the Heat Kernel preserves the constraint p < 1. As a 
consequence, we consider the modihed Fokker-Planck type equation 


(1-4) dtpt - Ap t + V • (pt.P a dm(pt) M) = 0, 

which can also be written equivalently for the variables (p,p) as 

n r x f d t pt - Apt + V • (pt{u t - Vpt)) = 0, in [0, T] x Cl, 

( p > 0, p < 1, p(l — p) = 0, in [0, T\ x Cl. 

As usual, these equations are complemented by no-hux boundary conditions and by an initial datum 
p(0, x) = po{x). 

Roughly speaking, we can consider that this equation describes the law of a motion where each 
agent solves the stochastic differential equation 

dX t = ( ut(X t ) - Xp t (Xt)) d t + y/2dB t . 


This last statement is just formal and there are several issues defining an SODE like this: indeed, the 
pressure variable is also an unknown, and globally depends on the law pt of fQ. Hence, if we wanted 
to see this evolution as a superposition of individual motions, each agent should somehow predict 
the evolution of the pressure in order to solve his own equation. This reminds of some notions from 
the stochastic control formulation of Mean-Field Games, as introduced by J.-M. Lasry and P.-L. 
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Lions, even if here there are no strategic issues for the players. For MFG with density constraints, 
we refer to [ 8 , 33, 39]. 

However, in this paper we will not consider any microscopic or individual problem, but only 
study the parabolic PDE (1.5). 

1.3. Structure of the paper and main results. The main goal of the paper is to provide 
an existence result, with some extra estimates, for the Fokker-Planck equation (1.5) via time 
discretization, using the so-called splitting method (the two main ingredients of the equation, 
i.e. the advection with diffusion on one hand, and the density constraint on the other hand, are 
treated one after the other). In Section 2 we will collect some preliminary results, including what 
we need from optimal transport and from the previous works about density-constrained crowd 
motion, in particular on the projection operator onto the set JC. In Section 3 we will provide the 
existence result we aim at, by a splitting scheme and some entropy bounds; the solution will be 
a curve of measures in HC 2 ([0, T]; W^H)) (absolutely continuous curves with square-integrable 
speed). In Section 4 we will make use of BV estimates to justify that the solution we just built 
is also Lip([0, T]; Wi(H)) and satisfies a global BV bound ||/Ot||_By < C (provided po € BV): this 
requires to combine BV estimates on the Fokker-Planck equation (which are available depending 
on the regularity of the vector field u ) with BV estimates on the projection operator on 1C (which 
have been recently proven in [14]). Section 5 presents a short review of alternative approaches, 
all discretized in time, but based either on gradient-flow techniques (the JKO scheme, see [18]) 
or on different splitting methods. Finally, in the Appendix A we detail the BV estimates on the 
Fokker-Planck equation (without any density constraint) that we could find; this seems to be a 
delicate matter, interesting in itself, and we are not aware of the sharp assumptions on the vector 
field u to guarantee the BV estimate that we need. 

2. Preliminaries 

2.1. Basic definitions and general facts on optimal transport. Here we collect some tools 
from the theory of optimal transportation, Wasserstein spaces, its dynamical formulation, etc. 
which will be used later on. We set our problem either in a compact convex domain C 
with smooth boundary or in the d —dimensional flat torus Cl := T d (even if we will not adapt all 
our notations to the torus case). We refer to [41, 40] for more details. Given two probability 
measures p, v € V{Cl) and for p > 1 we define the usual Wasserstein metric by means of the 
Monge-Kantorovich optimal transportation problem 

W p (p : u) := inf ( [ \x - y\ p dy(x, y) : 7 eU(p,u)\ P , 

Ul2x0 J 

where n(/u, v) := {7 € V(Cl X Cl) : {tt x )^ = p, ('K y )#'i = u} and n x and n y denote the canonical 
projections from Cl x Cl onto Cl. This quantity happens to be a distance on V(Cl) which metrizes 
the weak-* convergence of probability measures; we denote by W p (fl) := (V(Cl), W p ). i.e. the space 
of probabilities on Cl endowed with this distance. 

Moreover, in the quadratic case p = 2 and under the assumption p <C C d (the d —dimensional 
Lebesgue measure on Cl) in the late 80’s Y. Brenier showed (see [ 6 , 7]) that actually the optimal 7 
in the above problem (the existence of which is obtained simply by the direct method of calculus 
of variations) is induced by a map, which is the gradient of a convex function, i.e. there exists 
S : Cl —^ Cl and i/j : Cl —>• R convex such that S = and 7 := (id ,S)#p. The function if; is 
obtained as ip(x) = \\x\ 2 — </?(x), where tp is the so-called Kantorovich potential for the transport 
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from fi to v, and is characterized as the solution of a dual problem that we will not develop here. 
In this way, the optimal transport map S can also be written as S{x) = x — V</?(x). Later in the 
90’s R. McCann (see [34]) introduced a notion of interpolation between probability measures: the 
curve fit '■= {{T — t)x + ty)# 7, for t € [0, T] (T > 0 is given), gives a constant speed geodesic in 
the Wasserstein space connecting fi$ := fi and /it := v. 

Based on this notion of interpolation in 2000 J.-D. Benamou and Y. Brenier used some ideas 
from fluid mechanics to give a dynamical formulation to the Monge-Kantorovich problem (see [5]). 
They showed that 

T p_i Wgifi, v) = inf {B p (E, fi) : d t fi + V • E = 0, fi 0 = fi, fi T = v} ■ 


Here B p is a functional defined on pairs ( E,n ), where Li is a d-dimensional vector measure on 
[0, T] xO and fi = (fit)t. is a Borel-measurable family of probability measures on Q. This functional 
is defined to be finite only if E = Et <S> df (i.e. it is induced by a measurable family of vector measures 
on fl: /[ oT]xfi€(t’ x )' d E(t,x) = / Q T d t f n £(t,x)- dE t (x) for all test functions £ € C°([0,T] xH;M d )) 
and in this case it is defined through 


B p (E,fi) 


[ [ - \v t \ p dfi t (x)dt, 

Jo Jn P 

Too, 


if E t = v t - fit 
otherwise. 


It is well-known that B p is jointly convex and l.s.c. w.r.t the weak-* convergence of measures (see 
Section 5.3.1 in [40]) and that, if dtfi + V-E = 0, then B P (E , fi) < Too implies that t fit is a curve 
in Ah7 p ([0, T]; Wp(H)) 1 . In particular it is a continuous curve and the initial and final conditions 
on fiQ and fir are well-defined. 

Coming back to curves in Wasserstein spaces, it is well known (see [2] or Section 5.3 in [40]) that 
for any distributional solution fit (being a continuous curve in Wp(H)) of the continuity equation 
dtfi + V ■ E = 0 with Et = vt ■ fit, we have the relations 

W\w v {t) <\\v t \\ L ^ t and W p (fit, fi s ) < J \fi'\ Wp {r) dr, 

where we denoted by \fi'\w p {j) the metric derivative w.r.t. W p of the curve fit (see for instance [3] 
for general notions about curves in metric spaces and their metric derivative). For curves fit that 
are geodesics in W P (H) we have the equality 

W p (fi 0 ,fii) = \fi , \ Wp {t)dt = J \\vtWifcdt. 

The last equality is in fact the Benamou-Brenier formula with the optimal velocity field vt being 
the density of the optimal E t w.r.t. the optimal fi t . This optimal velocity field vt can be computed 
as vt := ( S — id) o ( St ) -1 , where St := (1 — t)id + tS is the transport in McCann’s interpolation 
(we assume here that the initial measure fio is absolutely continuous, so that we can use transport 
maps instead of plans). This expression can be obtained if we consider that in this interpolation 
particles move with constant speed S( x) — x, but x represents here a Lagrangian coordinate, and 
not an Eulerian one: if we want to know the velocity at time f at a given point, we have to find 
out first the original position of the particle passing through that point at that time. 


^Here AC V ($,T]\ Wj>(fi)) denotes the class of absolutely continuous curves in Wp(fl) with metric derivative in 
L p . See the connection with the functional B p . 
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In the sequel we will also need the notion of entropy of a probability density, and for any 
probability measure g € V(fl) we define it as 

g(x) log g(x) dx, if g <C C d , 

I +oo, otherwise. 

We recall that this functional is l.s.c. and geodesically convex in W^fl). 

As we will be mainly working with absolutely continuous probability measures (w.r.t. Lebesgue), 
we often identify measures with their densities. 

2.2. Projection problems in Wasserstein spaces. Our analysis strongly relies on the projec¬ 
tion operator Tfc in the sense of W^. Here K, := {p E T’(fl) : p < 1} and 

Pk[/A ■= argmin pe/c ^W 2 2 (//,p). 

We recall (see [30, 38] and [14]) the main properties of the projection P/c operator. 

• As far as is compact, for any probability measure p, the minimizer in min pg £ ^lE 2 (p, p) 
exists and is unique, and the operator P ^ is continuous (it is even C 0,1 / 2 for the W -2 distance). 

• The projection P/c[p\ saturates the constraint p < 1 in the sense that for any p € V(£l) 
there exists a measurable set B C Q such that Pjc [p] = 1 b + h ac l B c , where p ac is the 
absolutely continuous part of p. 

• The projection is characterized in terms of a pressure field, in the sense that p = Pjc [p] if 
and only if there exists a Lipschitz function p > 0, with p( 1 — p) = 0, and such that the 
optimal transport map S from p to p is given by S’ := id — S7p = id + 'S/p. 

• There is (as proven in [14]) a quantified BV estimate: if p G BV (in the sense that it is 
absolutely continuous and that its density belongs to BV(Sl)), then Pk\p\ is also BV and 

TV(p K [p],n)<TV(p,n)- 

This last BV estimate will be crucial in Section 4, and it is important to have it in this very 
form (other estimates of the form TV(Pjc[p], H) < aTV(p,Q) + b would not be as useful as this 
one, as they cannot be easily iterated). 

3. Existence via a splitting-up type algorithm (Main Scheme) 

Similarly to the approach in [31] (see the algorithm (13) and Theorem 3.5) for a general, non¬ 
gradient, vector field, we will build a theoretical algorithm, after time-discretization, to produce 
a solution of (1.5). Let us remark that splitting-type methods have been widely used in other 
contexts as well, see for instance the paper [10] which deals with splitting methods for Fokker-Planck 
equations and for more general gradient flows in metric and Wasserstein spaces, or [24] where a 
splitting-like approach is used to attack PDEs which are not gradient flows but “perturbations” of 
gradient flows. 

In this section the spontaneous velocity field is a general vector field u : [0, T] x Si —>• (not 
necessarily a gradient), which depends also on time. The only assumption we require on u is the 
following: 

(U) 



u G L°°([0, T] x fi; R d ). 
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We will work on a time interval [0, T] and in a bounded convex domain Q C W l (the case of the flat 
torus is even simpler and we will not discuss it in details). We consider po G 7 ?ac (fi) to be given, 
which represents the initial density of the population, and we suppose po G 1C. 


3.1. Splitting using the Fokker-Planck equation. Let us consider the following scheme. 


Main scheme: Let r > 0 be a small time step 
with N := [T/rJ. Let us set Pq := po and for 
every k G {1 ,,N} we define p k+ i from p T k in 
the following way. First we solve 
(3.1) 

f dtQt - A Qt + V ■ ( QtUt+kr) = o, t €]0,t], 

I A) = Pi, 

equipped with the no-flux boundary condition 
(pt(Vg t — ut) ■ n = 0 a.e. on dQ) and set 
Pk+ i = p ic[pl+ 1 ], where p T k+1 = g T . See Figure 
1 on the right. 



Figure 1. One time step 


Let us remark first that by classical results on parabolic equations (see for instance [25]), since 
u satisfies the assumption (U), Problem 3.1 admits a unique distributional solution. 

The above algorithm means the following: first follow the Fokker-Planck equation, ignoring the 
density constraint, for a time r, then project. In order to state and prove the convergence of the 
scheme, we need to define some suitable interpolations of the discrete sequence of densities that we 
have just introduced. 

First interpolation. We define the following curves of densities, velocities and momenta con¬ 
structed with the help of the p/Cs. First set 


/ 92 (t-fcr), if t € [kr, (k + 1 / 2 )r[, 

1 (id + 2 ((fc + 1)t — t)Vp^ +1 )_^ p k+l , if t G [(k + l/2)r, (A: + l)r[, 


where Qt is the solution of the Fokker-Planck equation (3.1) with initial datum p T k and VpA 1 arises 
from the projection of p T k+ \- more precisely (id + T'S7p T k+l ) is the optimal transport from p k+1 to 
Pk+i- What are we doing? We are fitting into a time interval of length r the two steps of our 
algorithm. First we follow the FP equation (3.1) at double speed, then we interpolate between the 
measure we reached and its projection following the geodesic between them. This geodesic is easily 
described as an image measure of p k+1 through McCann’s interpolation. By the construction it is 
clear that p/ is a continuous curve in V(iT) for t G [0, T], We now define a family of time-dependent 
vector fields though 


~ 2 ^ 2 q-fc!) ) + 2ut ’ if t G [kr, (k + l/2)r[, 

“ 2 V Pfc +1 ° ( id + 2 (( k + f) r “ f)Vp;[; +1 ) -1 , if t G [(k + 1/2)r, (k + l)r[, 


and, finally, we simply define the curve of momenta as EJ := p/v/. 
Second interpolation. We define another interpolation as follows. Set 


Pt '■= Qt-kr, if t G [kr, (k + l)r[, 
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where Qt is (again) the solution of the Fokker-Planck equation (3.1) on the time interval [0 ,t] with 
initial datum p T k . Here we do not double its speed. We define the curve of velocities 

~T ^ if t € [hr, ( k + 1 )r[, 

Qt—kr 

and we build the curve of momenta by EJ := pjvj ■ 

Third interpolation. For each r, we also define piecewise constant curves, 

Pt'-=Pk+n if t € [fcr, (A: + l)r[, 

v T t ■= Vpfc +1 , if i € [fcr, (A: + l)r[, 

and EJ := We remark that p£ +1 (l — p^ +1 ) = 0, hence the curve of momenta is just 

EJ '■= v Pfc+D if t € [&T, (A; + l)r[. 

Mind the differences in the construction of pj, pj and pj (hence in the construction of vj, vj 
and vj and EJ. EJ and EJ): 

1) the first one is continuous in time for the weak-* convergence, while the second and third ones 
are not; 

2) in the first construction we have taken into account the projection operator explicitly, while 
in the second one we see it just in an indirect manner (via the ‘jumps’ occurring at every time of 
the form t = kr). The third interpolation is piece-wise constant, and at every time it satisfies the 
density constraint; 

3) in the first interpolation the pair ( p T ,E T ) solves the continuity equation, while in the other 
two they do not. This is not astonishing, as the continuity equation characterizes continuous curves 
in W 2 (n). 

In order to prove the convergence of the scheme above, we will obtain uniform AC 2 ([0, T]; W 2 (O)) 
bounds for the curves p T . A key observation here is that the metric derivative (w.r.t. Wf) of the 
solution of the Fokker-Planck equation is comparable with the time differential of the entropy 
functional along the same solution (see Lemma 3.2). Now we state the main theorem of this 
section. 

Theorem 3.1. Let po € K, and u be a given desired velocity field satisfying (U). Let us consider the 
interpolations introduced above. Then there exists a continuous curve [0, T] 3 t i-a p t G W 2 (fi) and 
some vector measures E,E,E € 9Jt([0,T] x H) such that the curves p T , p T , p T converge uniformly 
in W 2 (H) to p and 

E T A E, E T ^ E, E T A E, in 9Jl([0, T\ x Ll) d , as r -> 0. 

Moreover E = E — E and for a.e. t € [0, T] there exist time-dependent measurable vector fields 
Vf,Vf, vt such that 

(1) E = pv, E = pv, E = pv, 

(2) [ (||vt |||2 + ||vt ||^2 + ll^illl 2 ) dt <+ 00 , 

J Q V Pt Pt Pt/ 

(3) v t = v t - v t pt - a.e., E t = p t u t - Vp t and v t = Vp £ , p t - a.e., 
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where p € L 2 ([ 0, T\; p > 0 and p{ 1 — p) = 0 a.e. in [0, T] x Q. As a consequence, the pair 

(p,p) is a weak solution of the problem 

d t pt - Ap t + V • (ptOt - Vpt)) = 0, in [0,T] x Q, 

( o Pt> 0, Pt < 1, Pt( 1 - Pt) = 0, in [0, T] x n, 

1 ' Ptiypt -u t + Vp t ) • n = 0, on [0, T] x <9fi, 

k P( 0 , ■) = Po- 

To prove this theorem we will use the following tools. 


Lemma 3.2. Let us consider a solution gt of the Fokker-Planck equation on [0,T] x fi with the 
velocity field u satisfying (U) and with no-flux boundary conditions on [0, T] x dLl. Then for any 
time interval ]a, 6 [ we have the following estimate 


(3.3) 



Vet 

Qt 


+ u t . 


Qt dx d t < £(g a ) - £{Qb) + 


1 


la 



ut\ 2 Qt dx d t 


In particidar this implies 


(3.4) 


\ [ \6t\w 2 dt ^ £ (ea)-£(0b) + \ [ [ \u t \ 2 Qtdxdt, 
1 Ja Z Ja Jn 


where |p(|yp 2 denotes the metric derivative of the curve t eA Qt € >V 2 (^). 


Proof. To prove this inequality, we will first make computations in the case where both u and q 
are smooth, and q is boundeded from below by a positive constant. In this case we can write 
d 


d t 


£(gt)= / (log Q t + l)<9tptdx = / log pt(Apt — V • (QtUt)) dx 

\VQt\ 2 


Qt 


+ u t ■ V Q t dx 


where we used the conservation of mass (i.e. J n dtQt dx = 0 ) and the boundary conditions in the 
integration by parts. We now compare this with 

2 , ^ /• /-, 1^7 i 2 


VQt 

Qt 


+ u t 


Qtdx - ^ / \u t \ 2 Qt dx = 


1 I VotY 

2 Qt 


- VQt ■ Ut dx 


< 


f 

In V Qt 


- Vpt - ut) dx = T(pt). 


This provides the first part of the statement, i.e. (3.3). If we combine this with the fact that the 
metric derivative of the curve 1 1 —)• Qt is always less or equal than the L 2 t norm of the velocity field 
in the continuity equation, we also get 


T:\Qt\w2 ~ o 



u t \ 2 Qt < 



and hence (3.4). 

In order to prove the same estimates without artificial smoothness and lower bound assumptions, 
we can act by approximation. We approximate the density Q a by smooth and strictly positive 
densities (by convolution, so that we guarantee in particular £(pk) —>• £(p a )), and the vector 
field u with smooth vector fields u k (strongly in L 4 ([a, 6 ] x fl), keeping the L°° bound). If we 
call Q k the corresponding solution of the Fokker Planck equation, it satisfies (3.3). This implies 
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a uniform bound (w.r.t. k) for y/ffi in L 2 ([a,b\] and hence a uniform bound on g k in 

L 2 ([a, 6 ] x Q). From these bounds and the uniqueness of the solution of the Fokker-Planck equation 
with L°° drift we deduce g k —>• g. The semicontinuity of the left-hand side in (3.3) and of the 
entropy term at t = b, together with the convergence of the entropy at t = a and the convergence 
la In \ uk \ 2 Q k dx d t — > f n \u\ 2 gdx dt (because we have a product of weak and strong convergence 
in L 2 ) allow to pass (3.3) to the limit. □ 

Corollary 3.3. From the inequality (3.4) we deduce that 

£(g b ) - £{g a ) < \ f [ \u t \ 2 g t dxdt, 
hence in particular for u satisfying (U) we have 

£(eb) -£{Qa) < (b — a). 

As a consequence, if g a < 1, then we have 

£(eb) < ^IMI 

The same estimate can he applied to the curve p T , with a = kr and b &}kr, (k + l)r[, thus obtaining 
£(pt) < Ct for every t. 


Lemma 3.4. For any p S V(£l) we have £ (Pjc[p]) ^ £{p)- 

Proof. We can assume p <C C d , otherwise the claim is straightforward. As we pointed out in Section 
2.2, we know that there exists a measurable set B C such that 


Pjc[p] = 1 B + P^-B c - 


Hence it is enough to prove that 


p log pdx >0 = / Pic[p] log P/c[p] dx, 


IB 


IB 


as the entropies on B c coincide. As the mass of p and P/c[p\ is the same on the whole and they 

coincide on B c , we have / p(x) dx= P)c[p\ dx = \B\. 

J B Jb 


Then, by Jensen’s inequality we have 
1 


LB 


IB 


1 


—— / p log p dx > I —— / pdx log —— / pdx = 0 . 


B 


IB 


1 


LB 


The entropy decay follows. 


□ 


To analyze the pressure field we will need the following result. 

Lemma 3.5. Let {p T } r >o a bounded sequence in L 2 ([0, T]; i/ 1 (H)) and {p T } r >o a sequence of 
piecewise constant curves valued in W 2 (H), which satisfy H 2 (p r (a), p T {b)) < Cy/b — a + r for all 
a < b E [0, T] for a fixed constant C. Suppose that 

P r > 0, p T { 1 - p T ) = 0, p T < 1, 

and that 

p T p weakly in B 2 ([0, T]; H l (H)) and p T —>• p uniformly in W^H). 

Then p{ 1 — p) = 0 a.e. in [0, T] x H. 
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The proof of this result is the same as in Step 3 of Section 3.2 of [30] (see also [37] and Lemma 
4.6 in [15]). We omit it in order not to overburden the paper. 

The reader can note the strong connection with the classical Aubin-Lions lemma [4], applied to 
the compact injection of L 2 into Lf . Indeed, from the weak convergence of p T in L\[0,T\-H l (Q)), 
we just need to provide strong convergence of p T in L 2 ([0, T]; H 1 (fl)). If instead of the quasi-Holder 
assumption of the above lemma we suppose a uniform bound of {p T } T in AC 2 {[ 0, T]; W 2 (H)) (which 
is not so different), then the statement can be really deduced from the Aubin-Lions lemma. Indeed, 
the sequence {p T } is bounded in L°°([0,T]-, L 2 (£l)) and its time-derivative would be bounded in 
L 2 ([0, T]; i7 _1 (H)). This strongly depends on the fact that the H~ l distance can be controlled 
by the W 2 distance as soon as the measures have uniformly bounded densities (see [29, 30]), a 
tool which also crucial in the proofs in [30, 37, 15]. Then, the Aubin-Lions lemma guarantees 
compactness in C°([0,T]; iL _1 (H)), which is more than what we need. 

Lemma 3.6. Let us consider the previously defined interpolations. Then we have the following 
facts. 

(i) For every t > 0 and k we have 

max { W 2 (Pk , Pk+i ), (PkiPk+i)} < tC ( £ (Pk) ~£(Pk+i)) +Ct 2 , 

where C > 0 only depends on ||rt|| - 

(ii) There exists a constant C, only depending on p 0 and HwIIl 00 ; such that B 2 (E T ,p T ) < C, 
B 2 {E T ,p T ) < C andB 2 (E T ,p T ) < C. 

(iii) For the curve [0, T] 3 t t-A pi we have that 

[ T \(PtY\w a *t<c, 

J 0 

for a C > 0 independent of r. Here we denoted by \{Pt)'\w 2 the metric derivative of the curve 
p T at t in W 2 . In particular, we have a uniform Holder bound on p T : W 2 (p T (a), p T (b)) < 
Cy/b — a for every b > a. 

(iv) E T ,E T ,E T are uniformly bounded sequences in 3Jl([0,T] x Q) d . 

Proof, (i) First by the triangle inequality and by the fact that p T k+1 = Pic[Pk+i\ we have that 

(3-5) W 2 {PkiPk+i) < W 2 (p k ,p k+ 1 ) + W 2 (pl +1 , p k+1 ) < 2W 2 (pl, pl + i)- 

We use (as before) the notation g t , t £ [0,r] for the solution of the Fokker-Planck equation (3.1) 
with initial datum p T k , in particular we have g T = p T k+l - Using Lemma 3.2 and since £o = P T k and 

Qr = Pfc i i we have by (3.4) and using W 2 (p k ,p k+1 ) < / \g[\w 2 dt 

Jo 

WiipbPk+i) < M (^J o \Pt\w 2 d^ J < 2 t(£:( £ > 0 ) ~£{g T )) + T J o J n \u k r+t\ 2 0 tdxdt 

< 2r ( £{p k ) — £{p k+ 1 )) + Ct 2 < 2t ( £{p k ) — £(pl+1 )) + Cr 2 , 

where C > 0 is a constant depending just on ||u||l°o. We also used the fact that £(p k+1 ) < £(p k+1 ), 
a consequence of Lemma 3.4. 

Now by the means of (3.5) we obtain 

Wi(p T k ,p T k+1 ) < tC(£(pD-£(pI + 1 ))+Ct 2 . 
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(ii) We use Lemma 3.2 on the intervals of type [kr, (k + l/2)r[ and the fact that on each interval 
of type [(k + l/2)r, (k + l)r[ the curve pf is a constant speed geodesic. In particular, on these 
intervals we have 

\{p t )'\w 2 = HWl\ = 2t \\VpI+iWl\ = 2W 2(pI+i,PI + i)- 

p t p l+1 

On the other hand we also have 

' 2 ||VpI +1 ||k = Wi(pl +1 ,~ P T k+ 1 ) < Wi(pl,fT k+ 1 ) < rC{£(pl)-£(pl +1 )) +Ct 2 . 


T 


H k+1 


Hence we obtain 

»(fc+l)T 


J kr 


Kll h( P i) dt 


/•(fe+l/2)r r 

J kr Jfl 


^ @2(t—kr) 
Q2(t—kr) 




Q2{t-kT) i x )dx dt + 4 


r(k-\-V)T 

' (k+l/2)r J 


|Vp]fe +1 | p T k+1 dxdt 


H k+1 


< C (£{pl) ~ £(pl+ 1)) + Ct + 2r||Vpfc +1 ||^2 

p l 

<C{£{pl)-£{pl + 1 ))+Cr. 

Hence by adding up we obtain 

B 2 (E\p T ) <Y,{C {S(p T k ) ~ £(pl+1 )) + Cr} = C {£{pl) - £(p T N+1 )) + CT<C. 


The estimate on B 2 (E T ,p T ) and B 2 (E T , p T ) are completely analogous and descend from the 
previous computations. 

r T 

(iii) The estimate on B 2 (E T , p T ) implies a bound on / | {pY)'\w 2 because v T is a velocity field 

Jo 

for p T (i.e., the pair (E T , p T ) solves the continuity equation). 

(iv) In order to estimate the total mass of E we write 


|£ r |([0,T] x ft) = [ [ \v T t \p T t dxdt 

Jo Jn 


< 


\vt\ 2 Pt dx 


pl dx dt 


< 




vf\ pfdxdtj < C. 

n ) 


The bounds on E T and E T rely on the same argument. 


□ 


Proof of Theorem 3.1. We use the tools from Lemma 3.6. 

Step 1. By the bounds on the metric derivative of the curves pf we get compactness, i.e. there 
exists a curve [0, T] 3 t eA pt E V(£l) such that p T (up to subsequences) converges uniformly in 
[0, T] w.r.t. W 2 , in particular weakly-* in V(£l) for all t E [0, T]. It is easy to see that p T and p T 
are converging to the same curve. Indeed we have pj = p T ~^ and p[ = p)W for |s(t) — t\ < r and 

| s(t) — t\ < r, which implies W 2 (pf , pf),W 2 (pJ ", pj) < Ct T This provides the convergence to the 
same limit. 

Step 2. By the boundedness of E T ,E T and E T in 9Jt([0, T] x £l) d we have the existence of 
E,E,E E 9Jt([0, T] x £l) d such that (up to a subsequence) E T -E E, E T e, E t E as r —* 0. 
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Now we show that E = E — E. Indeed, let us show that for any test function / € Lip([0, T] x Q) d 
we have 

rT 


fi 


f t .(El-(El + EJ))( dx, dt) 


0 , 


as r —> 0. First for each k € {0,..., N} we have that 


/■(fe+l/2)r r r(K+L)T r 

/ / ft ■ dx, dt) = / / f(t+kr)/2 ■ (-Vft_jfc T + «tet-fcr)(dx, dt) 

«/ kr J J kr J fl 


c(fe+l)r 


' fcr 

/■(fe+l)r /• _ r(K+t)T r 

/ ft-Eli dx, dt) + / / (f{t+kr)/2 - ft) ■ Ef (dx, dt) 

J kr J fl J kr J £7 


r(fc+l)r 


and 

r-(fc+l)r 


/•(K-tl)T r r(k+t)T r 

/ ffEJ{dx,dt) = / -f(t+(k+i)T)/ 2 °(. [ d + ((k + l)T-t)X7pl +1 )-X7p T k+1 pl +1 (dx,dt) 

J(k+l/2)rJn J kr Jfl 

ft ■ Ef (dx, dt) 

MK+ijT r 

+ (ft- f(t+(k+ i)r)/2 ° (id + ((A; + l)r - f))) • D t r p[ (dx, dt) 

J kr J fl 


' kr J 

r(k+l)T 

f kr 

r(k+ i)t 


This implies that 

rT 


[ [ ff(E x 
Jo Jn 


J-EZ + EZ)( dx, dt) 


< 


E 


-»(fc+l)7 


kr 


Lip(/) r / lA T |(dx, dt) 


/n 


/■(fe+l)T r 

+ E / Li p(/) r / (! + |w[|)|-Ei T |(dx, dt) 

7 J kr J Q, 


<TCUp(f)(\E T \(l0,T}xn) + \E T \([0,T}xn) + B 2 (E,p 
< rCLip(/), 


for a uniform constant C > 0. Letting r->0we prove the claim. 

.Step 5. The bounds on B 2 (E T , p T ), B 2 (E T , p T ) and B 2 (E T , p T ) pass to the limit by semicontinuity 
and allow to conclude that E, E and E are vector valued Radon measures absolutely continuous 
w.r.t. p. Hence there exist vt,vt,vt such that E = pv, E = pv and E = pit. 

Step 4■ We now look at the equations satisfied by E, E and E. First we use dtp T + V • E T = 0, 
we pass to the limit as r —>• 0, and we get 


d t p + V • E = 0. 

Then, we use E T = —Vp r + utp T , we pass to the limit again as r —>• 0, and we get 

E = —Vp + utp. 


To justify the above limit, the only delicate point is passing to the limit the term utp T , since u 
is only L°°, and p T converges weakly as measures, and we are a priori only allowed to multiply it 
by continuous functions. Yet, we remark that by Corollary 3.3 we have that £{pt) < Ct for all 
t £ [0, T\. In particular, this provides, for each t, uniform integrability for pf and turns the weak 
convergence as measures into weak convergence in L 1 . This allows to multiply by ut in the weak 
limit. 
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Finally, we look at E T . There exists a piecewise constant (in time) function p T (defined as p £ +1 
on every interval ]kr, (k + l)r]) such that p T > 0 , p T { 1 — p T ) = 0 , 

(3.6) ( f |Vp r | 2 ( dx, dt) = f f \\7p T \ 2 p T (dx, dt) = f f \v T \ 2 p T (dx, dt) < C 

Jo Jn Jo Jn Jo Jn 

and E T = X7p T p T = Vp T . The bound (3.6) implies that p T is uniformly bounded in L 2 (0, T; H 1 (0)). 
Since for every t we have |{pi = 0}| > \{p[ < 1} | > |I7| — 1, we can use a suitable version of 
Poincare’s inequality, and get a uniform bound in L 2 ([0, T\; L 2 (I7)) = L 2 ([0, T] x Q). Hence there 
exists p E L 2 ([0, T] x 17) such that p T —^ p weakly in L 2 as t —^ 0. In particular we have E = Vp. 
Moreover it is clear that p > 0 and by Lemma 3.5 we obtain p(l — p) = 0 a.e. as well. Indeed, the 
assumptions of the Lemma are easily checked: we only need to estimate W 2 (p T (a), p T (b)) for b > a, 
but we have 

W 2 (p T (a),p T (b)) = W 2 (p T (k a r), p T {k b r)) < Cy/k b - k a , for k b r < b + r and k a > a. 

Once we have E = Vp with p(l — p) = 0, p E L 2 ([0, T]\ LI 1 (17)) and p E L°°, we can also write 

E = Vp = pVp. 

If we sum up our results, using E = E — E, we have 

dtp — A p + V • (p(u — Vp)) = 0 together with p > 0, p < 1, p(l — p) = 0 a.e. in [0, T] x 17. 

As usual, this equation is satisfied in a weak sense, with no-flux boundary conditions. □ 

4. Uniform Lip([0,T];Wi) and BV estimates 

In this section we provide uniform estimates for the curves p T ,p T and p T of the following form: 
we prove uniform BV (in space) bounds on p T (which implies the same bound for p T ) and uniform 
Lipschitz bounds in time for the W\ distance on p T . This means a small improvement compared 
to the previous section in what concerns time regularity, as we have Lipschitz instead of AC 2 , 
even if we need to replace W 2 with W\. It is also important in what concerns space regularity. 
Indeed, from Lemma 3.2 one could deduce that the solution p of the FP equation (1.5) satifies 
y/p E L 2 ([0, T]; L7 1 (17)) and, using p < 1, also p E L 2 ([0, T]; 7L 1 (17)). Yet, this is just an integrable 
estimate in t , while the BV estimate of this section is uniform in the time variable. 

Nevertheless there is a price to pay for this improvement: we have to assume higher regularity 
for the velocity field. These uniform-in-time VFi-Lipschitz bounds are based both on BV estimates 
for the Fokker-Planck equation (see Lemma A.l from Appendix A) and for the projection operator 
Pic (see [14]). The assumption on u is essentially the following: we need to control the growth of 
the total variation of the solutions of the Fokker-Planck equation (3.1), and we need to iterate this 
bound along time steps. 

We will discuss in the Appendix the different BV estimates on the Fokker-Planck equation that 
we were able to find. The desired estimate is true whenever [['Ut||cM(n) is uniformly bounded 
ut ■ n = 0 on <917. It seems to be an open problem to obtain similar estimate under the only 
assumption that u is Lipschitz continuous. Of course, we will also assume po E BV(£l). Despite 
these extra regularity assumptions, we think these estimates have their own interest, exploiting some 
finer properties of the solutions of the Fokker-Planck equation and of the Wasserstein projection 
operator. 

Before entering into the details of the estimates, we want to discuss why we concentrate on BV 
estimates (instead of Sobolev ones) and on W\ (instead of W p , p > 1). The main reason is the role 
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of the projection operator: indeed, even if p E kF 1,p (12), we do not have in general Pjc\p] E W 1,p 
because the projection creates some jumps at the boundary of {Pk[p\ = 1}- This prevents from 
obtaining any W l,p estimate for p > 1. On the other hand, [14] exactly proves a BV estimate on 
P/C [p] and paves the way to BV bounds for our equation. Concerning the regularity in time, we 
observe that the velocity field in the Fokker-Planck equation contains a term in ^7 p/p. Since the 
metric derivative in W p is given by the L p norm (w.r.t. pt) of the velocity field, it is clear that 
estimates in W p for p > 1 would require spatial W 1)P estimates on the solution itself, which are 
impossible for p > 1 in this splitting scheme. We underline that this does not mean that uniform 
W 1,p are impossible for the solution of (1.5); it only means that they are not uniform along the 
approximation that we used in our Main Scheme to build such a solution. 

The precise result that we prove is the following. 

Theorem 4.1. Let us suppose that ||c'i-i — C an d Po E BV(Q). Then using the notations from 
the Main scheme and Theorem 3.1 one has ||p[||si/ < C and W\(p k , p k+1 ) < Ct. As a consequence 
we also have p E Lip([0, T]; Wi) D L°°([0, T]; BV (0)). 

To prove this theorem we need the following lemmas. 

Lemma 4.2. Suppose 11it*11Lip < C and ut ■ n = 0 on <912. Then for the solution g of (A.l) with 

velocity field v = u we have the estimate 

ii ii ii ii Ct 

||0t||L°° < ||0o||L°°e , 

where C = ||V • ut ||l°°■ 

Proof. Standard comparison theorems for parabolic equations allow to prove the results once we 
notice that f(t,x ) := HftolU 00 ^* is a supersolution of the Fokker-Planck equation, i.e. 

dtft > A/ t - V ■ {f t ut). 

Indeed, in the above equation the Laplacian term vanishes as / is constant in x , dtft = Cft and 
V • (ftUt) = ,/#V • ut + V/t • ut = /iV • ut < Cft where C = ||V • ut\\L°°- From this inequality, and 
from po < /o, we deduce pt < ft for all t. □ 

We remark that the above lemma implies in particular that after every step in the Main scheme 
we have p T k+1 < e TC < 1 + Ct, where c := ||V • u\\l°°- Let us now present the following lemma as 
well. 

Corollary 4.3. Along the iterations of our Main scheme, for every k we have Wi(p T k+l , p T k+l ) < tC 
for a constant C > 0 independent of r. 

Proof. With the saturation property of the projection (see Section 2.2 or [14]), we know that there 
exists a measurable set B C 11 such that p T k+l = p k+ itB + 1ms- O n the other hand we know that 

w i(Pk+nPk+ 1) = SU P / f(Pk+i -Pl+i)dx 

/GLipi (^), 0</<diam(r2) JQ 

= sup / f(Pk +1 — 1) dx < tC |n|diam(n). 

fEhipi (Q), 0</<diam(Q) J£l\B 

We used the fact that the competitors / in the dual formula can be taken positive and bounded 
by the diameter of 12, just by adding a suitable constant. This implies as well that C is depending 
on c, |0| and diam(12). □ 
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Proof of Theorem f.l. First we take care of the BV estimate. Lemma A.l in the Appendix guar¬ 
antees, for t g]/ct, (/c + 1)t[, that we have TV (pf) < CT+e CT TV{p T k ). Together with the BV bound 
on the projection that we presented in Section 2.2 (taken from [14]), this can be iterated, providing 
a uniform bound (depending on TV (po)- T and sup f \\ut Hcm) on ||p^||sv- Passing this estimate to 
the limit as r —> 0 we get p € L°°([0, T]; BV(VL)). 

Then we estimate the behavior of the interpolation curve p T in terms of W \. We estimate 

r{k+l)r f(k+ l)r r /ly^rl \ 

W\{p T k;pl+\) < I mr\ Wl dt< / / l—=^ + \ut\)ptdxdt 

Jkr Jkr JQ \ Pt J 

/»(/e+l)r 

< / \\pI\\bv dt + Ct < Ct. 

J kr 

Hence, we obtain 

w i(fibPk+ 1) < ^i(pLPfc+i) + tWfc+DPfc+i) < tC - 

This in particular means, for b > a, 

Wi(p T (a), p T (b)) < C(b -a + r). 

We can pass this relation to the limit, using that, for every t. we have p[ —>• pt in (and hence 

also in Wi(fi), since W\ < W 2 ), we get 

Wi{p(a),p(b)) < C(b-a ), 

which means that p is Lipschitz continuous in Wi (O). □ 

5. Variations on a theme: some reformulations of the Main scheme 

In this section we propose some alternative approaches to study the problem (1.5). The general 
idea is to discretize in time, and give a way to produce a measure p T k+l starting from p k . Observe 
that the interpolations that we proposed in the previous sections p T ,p T and p T are only technical 
tools to state and prove a convergence result, and the most important point is exactly the definition 
of Pk+V 

The alternative approaches proposed here explore different ideas, more difficult to implement 
than the one that we presented in Section 3, and/or restricted to some particular cases (for instance 
when u is a gradient). They have their own modeling interest and this is the main reason justifying 
their sketchy presentation. 

5.1. Variant 1: transport, diffusion then projection. We recall that the original splitting 
approach for the equation without diffusion ([31, 37]) exhibited an important difference compared 
to what we did in Section 3. Indeed, in the first phase of each time step (i.e. before the projection) 
the particles follow the vector field u and p k+ i was not defined as the solution of a continuity 
equation with advection velocity given by ut, but as the image of p k via a straight-line transport: 
p k+1 '■= (id + TUkr)#Pk■ One can wonder whether it is possible to follow a similar approach here. 

A possible way to proceed is the following: take a random variable X distributed according to 
p k , and define p k+ i as the law of X + ruk T (X) + B r , where B is a Brownian motion, independent 
of X. This exactly means that every particle moves starting from its initial position X, following 
a displacement ruled by u, but adding a stochastic effect in the form of the value at time r of a 
Brownian motion. We can check that this means 

Pk +1 : = Vr * ((id + TU kr )#Pk ) , 
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where r] T is a Gaussian kernel with zero-mean and variance r, i.e. r) T {x) := 
Then we define 


1 _i£p 

- c 4r 

(4 r7r )d/2 


Pk+1 '= P IC [pk+ 1 ] • 

Despite the fact that this scheme is very natural and essentially not that different from the 
Main scheme , we have to be careful with the analysis. First we have to quantify somehow the 
distance W p (p k , p k+1 ) for some p > 1 and show that this is of order r in some sense. Second, we 
need to be careful when performing the convolution with the heat kernel (or adding the Brownian 
motion, which is the same): this requires either to work in the whole space (which was not our 
framework) or in a periodic setting (D = T d , the flat torus, which is qutie restrictive). Otherwise, 
the “explicit” convolution step should be replaced with some other construction, such as following 
the Heat equation (with Neumann boundary conditions) for a time r. But this brings back to a 
situation very similar to the Main scheme , with the additional difficulty that we do not really have 
estimates on (id + TUk T )#P T k - 


5.2. Variant 2: gradient flow techniques for gradient velocity fields. In this section we 
assume that the velocity field of the population is given by the opposite of the gradient of a function, 
ut = —VVt a typical example is given when we take for V the distance function to the exit (see the 
discussions in [30] about this type of question). We start from the case where V does not depend 
on time, and we suppose V € VF 1,1 (D). In this particular case - beside the splitting approach - 
the problem has a variational structure, hence it is possible to show the existence by the means of 
gradient flows in Wasserstein spaces. 

Since the celebrated paper of Jordan, Kinderlehrer and Otto ([18]) we know that the solutions 
of the Fokker-Planck equation (with a gradient vector field) can be obtained with the help of 
the gradient flow of a perturbed entropy functional with respect to the Wasserstein distance W 2 ■ 
This formulation of the JKO scheme was also used in [30] for the first order model with density 
constraints. It is easy to combine the JKO scheme with density constraints to study the second 
order/diffusive model. As a slight modification of the model from [30], we can consider the following 
discrete implicit Euler (or JKO) scheme. As usual, we fix a time step r > 0, pj = po and for all 
k E { 1,2, , [N/t\} we just need to define p k+1 - We take 

(5.1) p T k+l = argmin pe7 , (n) {Y V ( X M X ) dx + s (p) + Mp) + -^W 2 2 (p>pD} , 

where I/c is the indicator function of /C, which is 


Ik{x) 


0, if x E /C, 
+ 00 , otherwise. 


The usual techniques from [18, 30] can be used to identify that System (1.5) is the gradient flow 
of the functional p i-a J(p) := / V(x)p(x) dx + £{p) + Ifc(p) and that the above discrete scheme 

Jn 

converges (up to a subsequence) to a solution of (1.5), thus proving existence. The key estimate 
for compactness is 

—W 2 2 (p^ + i,Pfc) < J(Pk) - J(Pk+ 1)5 

which can be summed up (as on the r.h.s. we have a telescopic series), thus obtaining the same 
bounds on B 2 that we used in Section 3. 
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Note that whenever D 2 V > XI, the functional p i-a f Q V ( x)p(x ) d x+£(p)+Ijc{p) is A-geodesically 
convex. This allows to use the theory in [2] to prove not only existence, but also uniqueness for 
this equation, and even stability (contractivity or exponential growth on the distance between two 
solutions) in W 2 - Yet, we underline that the techniques of [16] also give the same result. Indeed, 
[16] contains two parts. In the first part, the equation with density constaints for a given velocity 
field u is studied, under the assumption that — u has some nronotonicity properties: (— ut(x) + 
u t (y)) ■ (x — y) > X\x — y\ 2 (which is the case for the gradients of A-convex functions). In this 
case standard Gronwall estimates on the W 2 distance between two solutions are proved, and it is 
not difficult to add diffusion to that result (as the Heat kernel is already contractant in W 2 ). In 
the second part, via different techniques (mainly using the adjoint equation, and proving somehow 
L 1 contractivity), the uniqueness result is provided for arbitrary L°° vector fields u, but with the 
crucial help of the diffusion term in the equation. 

It is also possible to study a variant where V depends on time. We assume for simplicity 
that V € Lip([0,T] x O) (this is a simplification; less regularity in space, such as W 1,1 , could be 
sufficient). In this case we define 

Jt(p) '■= [ V t (x)p(x)dx + £(p) + I K (p) 

Jn 

and 

(5-2) p T k+ 1 = argmin peP(n) j J kT {p) + ^-\V$(p, p)T) j , 

The analysis proceeds similarly, with the only exception that the we get 

—W 2 (pI +1 ,pI) < J kr {p T k ) — J kr {p k+1 ), 

which is no more a a telescopic series. Yet, we have J kr {p T k _ |_i) > J(k+i)r(Pk+i ) + Lip(V)r, and we 
can go on with a telescopic sum plus a remainder of the order of r. In the case where ut is the 
opposite of the gradient of a A-convex function Vt, one could consider approximation by functions 
which are piecewise constant in time and use the standard theory of gradient flows. 

Let us remark here that the recent paper [1] gives another approach to deal with first order 
crowd motion models as limit of nonlinear-diffusion equations with gradient drift. This approach 
could be plausible also in the case when we add a simple diffusion term in the models studied in 
[!]• 

5.3. Variant 3: transport then gradient flow-like step with the penalized entropy func¬ 
tional. We present now a different scheme, which combines some of the previous approaches. It 
could formally provide a solution of the same equation, but presents some extra difficulties. 

We define now p k+1 := (id + ru kT )#p k and with the help of this we define 

Pk+i ■= argmin peK S(p) + ^W^(p,pj. + l ). 

In the last optimization problem we minimize a strictly convex and l.s.c. functionals, and hence we 
have existence and uniqueness of the solution. The formal reason for this scheme being adapted to 
the equation is that we perform a step of a JKO scheme in the spirit of [18] (without the density 
constraint) or of [30] (without the entropy term). This should let a term —A p — V • (p'V'p) appear 
in the evolution equation. The term V • ( pu ) is due to the first step (the definition of p k+ i)- To 
explain a little bit more for the unexperienced reader, we consider the optimality conditions for the 
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above minimization problem. Following [30], we can say that p £ /C is optimal if and only if there 
exists a constant and a Kantorovich potential <p for the transport from p to pj, such that 


1 on (hip + < £, 

p = 0 on (in p + £) > l, 

^ £ [0, 1 ] on (in p+^)=£. 

We then define p = (£ — lnp — £) + and we get 
p £ press(p). Moreover, p— a.e. Vp = — ^T — 

We then use the fact that the optimal transport 
is of the form T = id—V<p and obtain a situation 
as is sketched in Figure 2. 


Pk' 


Pk+i 



id -r(u {fc+ i) T —Vp—^)+o(r) 


Figure 2. One time step 


'Pk +1 


Notice that (id + Tu kT ) 1 o (id + r(Vp + Vp/p)) = id — t(u^+i)t ~ V p — Vp/ p) + o(r) provided 
u is regular enough. Formally we can pass to the limit r —>• 0 and have 


d t p - Ap + V • (p(u - Vp)) = 0. 

Yet, this turns out to be quite naive, because we cannot get proper estimates on W 2 (p k , p/ +1 ). 
Indeed, this is mainly due to the hybrid nature of the scheme, i.e. a gradient flow for the diffusion 
and the projection part on one hand and a free transport on the other hand. The typical estimate 
in the JKO scheme comes from the fact that one can bound H^p^, p£ +1 ) 2 / r with the opposite of 
the increment of the energy, and that this gives rise to a telescopic sum. Yet, this is not the case 
whenever the base point for a new time step is not equal to the previous minimizer. Moreover, the 
main difficulty here is the fact that the energy we consider implicitly takes the value +oo, due to 
the constraint p £ /C, and hence no estimate is possible whenever p / +1 ^ 1C. As a possible way to 
overcome this difficulty, one could approximate the discontinuous functional Ijc with some finite 
energies of the same nature (for instance power-like entropies, even if the best choice would be an 
energy which is Lipschitz for the distance W 2 ) ■ These kinds of difficulties are matter of current 
study, in particular for mixed systems and/or multiple populations. 


Appendix A. BV- type estimates for the Fokker-Planck equation 

Here we present some Total Variation (TV) decay results (in time) for the solutions of the 
Fokker-Planck equation. Some are very easy, some trickier. The goal is to look at those estimates 
which can be easily iterated in time and combined with the decay of the TV via the projection 
operator, as we did in Section 4. 

Let us take a vector field v : [0,+oo[xH —» M. d (we will choose later which regularity we need) 
and consider in H the problem 

( d t p t - Ap t + V • ( p t v t ) = 0, in ]0, +oo[xH, 

(A.l) < p t (Vp t - v t ) ■ n = 0, on [0, +oo[x<9fl, 

l p( 0 ,-) = po, in Q, 

for po € bv(Q) nv(n). 

Lemma A.l. Suppose H^t lloi’i — C f or a H t £ [0,+oo[. Suppose that either H = T d or that Q is 
convex and v ■ n = 0 on dd. Then, we have the following total variation decay estimate 

(A. 2 ) f |Vpt| dx < C(t — s) + e c ^~ sS> f |Vp s |dx, V0 <s<t, 

Jo. Jfl 

where C > 0 is a constant depending just on the C 1,1 norm of v. 
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Proof. First we remark that by the regularity of v the quantity ||u|| L oo + \\Dv\\ L oo + || V(V • u)||^oo 
is uniformly bounded. Let us drop now the dependence on t in our notation and calculate in 
coordinates 


di./„ |V ' ,|dx = 


Vp 


|Vp| 

= -/E 


■ V(<9 t p) dx = 


Vp 

|Vp| 


■ v(Ap - v ■ (■ vp))dx = I Y -£h ( 5Z Kv ~ ( v ' ( w ^)) 


/n“ |Vp| 


Pj PkPkiPij 


dx + B 1 - I Y fph ( v ijP + V iPj + v jPi + v% Pij) dx 




|Vp| |Vp| 3 

<B 1 + C + C [ |Vp| da; + [ |Vp| |V ■ v\ dx + B 2 
J n Jn 

V B\ + i?2 + C + C f |Vp| dx. 

Jn 

Here the B, are the boundary terms, i.e. 

‘ Pjfd Pij 


n-| y p| 


By : = 


|Vp| 


d^ 


rl— 1 


and H 2 := — 


/( 


u • n)| Vp| d^ 


d—1 


The constant C > 0 only depends on ||u||i,oo + ||V • v \\ l °° + ||V(V ■ v )\\ l °°. We used as well the fact 

\( Pij _ Pj PkPkiPij ) j ^ n 
2-^ IVnl l\7 n|3 


that — 


^ pj.A; 


JVpI |Vp|= 

Now, it is clear that in the case of the torus the boundary terms B\ and B 2 do not exist, hence 
we conclude by Gronwall’s lemma. In the case of the convex domain we have B 2 = 0 (because of 
the assumption v ■ n = 0) and B\ < 0 because of the next Lemma A.2. □ 


Lemma A.2. Suppose that u : H —>• is a smooth vector field with u ■ n = 0 on dfl, p is a 
smooth function with Vp • n = 0 on dfl, and that H C is a smooth convex set that we write as 
fl = {h < 0} for a smooth convex function h with |V/i| = 1 on dfl (so that n = V/i on dfl). Then 
we have, on the whole boundary dfl, ’Y^ v ^jPj n% = ~ u l hjjPj. 

i,j i,j 

In particular, we have Y^PdPj nl < 0. 

hj 


Proof. The Neumann boundary assumption on u means 11 ( 7 (t)) ■ \/h( r y(t)) = 0 for every curve 7 
valued in dfl and for all t. Differentiating in t, we get 

+ X)«*(7(0) /l ii(7(0)(7 , ( i )) J = °- 

i,j i,j 

Take a point xq E dfl and choose a curve 7 with 7 (fo) = xq and 7 / (io) = Vp(x 0) (which is possible, 
since this vector is tangent to dfl by assumption). This gives the first part of the statement. The 
second part, i.e. Y^ PijPj ra * — 0, is obtained by taking u = Vp and using that D 2 h(x 0 ) is a positive 

i,j 

definite matrix. □ 


Remark A.3. If we look attentively at the proof of Lemma A.l, we can see that we did not really 
exploit the regularizing effects of the diffusion term in the equation. This means that the regularity 
estimate that we provide are the same that we would have without diffusion: in this case, the 
density pt is obtained from the initial density as the image through the flow of v. Thus, the density 


dx 
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depends on the determinant of the Jacobian of the flow, hence on the derivatives of v. It is normal 
that, if we want BV bounds on p t , we need assumptions on two derivatives of v. 


We would like to prove some form of BV estimates under weaker regularity assumptions on v, 
trying to exploit the diffusion effects. In particular, we would like to treat the case where v is only 
C 0,1 . As we will see in the following lemma, this degenerates in some sense. 


Lemma A.4. Suppose that Q is either the torus or a smooth convex set = {h < 0} parameterized 
as a level set of a smooth convex function h. Letvt : II —> W l be a vector field for t € [0, T], Lipschitz 
and bounded in space, uniformly in time. In the case of a convex domain, suppose v ■ n = 0 on 8Tl. 
Let H :W l —>• M be given by H(z) := \J e 2 + \z\ 2 . Now let p t (sufficiently smooth) be the solution of 
the Fokker-Planck equation with homogeneous Neumann boundary condition. 

Then there exists a constant C > 0 (depending on v and Q ) such that 


(A.3) f H(Vp t ) dx < [ H(Vp 0 )dx + C£t + — [ ||p s |||<x» ds. 

Jn Jn £ Jo 


Proof. First let us discuss about some properties of H. It is smooth, its gradient is VH(z) 
and it satisfies S7H(z) ■ z < H(z), \/z € M d . Moreover its Hessian matrix is given by 


z 

W) 




6 ij H 2 (z) - zV 
H 3 (z) 


Id ~ 


.</( H(z) H 3 (z) 


• z, Vze 


f 1 if i == i 

where 5 13 = | q' if j ^ j K ronec ^ er symbol. Note that, from this computation, the 

matrix D 2 H > 0 is bounded from above by —, and hence by e _1 . Moreover we introduce a 

H 

uniform constant C > 0 such that ||v||ic»|fi| + ||V • v\\l°° + ||Uu||loo < C. 

d /* 

Now to show the estimate of this lemma we calculate the quantity — / H(Vpt) da™. 


d 

dt 


H{N Pt ) dx= [ VH(Vp t )-d t Vpt dx= [ VH(Vp t ) ■ \7(Ap t — V • (v t pt)) dx 
! Jci Jci 

= [ \7H(\7p t )-\7Ap t dx- f VFT(Vpt) • V(V • (u t p t )) dx 
Jfl J Q 

=■ (I) + (II) 
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Now we study each term separately and for the simplicity we drop the t subscripts in the followings. 
We start from the case of the torus, where there is no boundary term in the integration by parts. 


(I) 



VH(Vp) • VApdx 


j 

Jn 


^ Hj(Vp)pju dx 





(II) 


- f VH(Vp) • V(V • (yp)) dx = - / V H 3 (V p)(v l p) VJ dx 
Jn Jn ij 

/ y Hjk (V p)Pkiv)p dx + / yH jk (Vp)p k ypjdx 

I" 7Tu TTi 


' i,j,k 

{I I a ) + (II b ). 


' i,j,k 


First look at the term ( II a ). Since the matrix Hj k is positive dehnite, we can apply a Young 
inequality for each index i and obtain 



< \m+C\\p\\ 2 L 4D 2 H\\ L ~. 


The L 2 norm in the second term will be estimated by the L°° norm for the sake of simplicity (see 
Remark A.5 below). 

For the term (lib) we first make a point-wise computation 


y H jk (\7p)p ki v l pj 
i,j,k 


H 3 (Vp) 4 


yiDfp • ( e 2 I d + \Vp\ 2 I d -X7p® Vp) • Vp]i 


H 3 (Vp) 


y v i D 2 i p-vp=-e 2 y v i d i 


H(y p) J 


where D 2 p denotes the i th row in the Hessian matrix of p and we used (| Vp\ 2 Id — Vp <8> Vp) -Vp = 0. 
Integrating by parts we obtain 

(lib) = e 2 J (V • dx < C£ 2 \\ 1 /H\\l°° < Ce, 

where we used H(z) > e. 

Summing up all the terms we get and using ||D 2 H|| < e _1 we get 


4 [ H(yPt)&x < — ^|(/)| + C\\pt\\y \\D 2 H\\l°° + Ce < Ce + C\\pt\\ 2 Lao e~ l , 
dt Jn z 

which proves the claim. 

If we switch to the case of a smooth bounded convex domain Q, we have to handle boundary 
terms. These terms are 


Idn 


y Hj(vp) Pij 


n — 




tan 


yHj(v P )p 


v j n 


hi 
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where we ignored those terms involving n * l v' (i.e., the integration by parts in (lib), and the term 
Hj(yp) Pj n l v l in the integration by parts of (II a )), since we already supposed v ■ n = 0. We use 
here Lemma A.2, which provides 


p)p ij n l -pH j {Wp)v)n l = R {pjPijri - ppjv)n l ) = 

ij ' i,j 


HiS p) 


^ ' {Pjhijpt ppjhijv) . 


hj 


If we use the fact that the matrix D 2 h is positive definite and a Young inequality, we get ■ PjhijPi > 
0 and 

P \Pj h v yl \ Pj h vPi + P 2 v j h ij v\ 

i,j i,j hi 


which implies 


1 

hNp ) 


J2 (P.iP'.i n ‘ 

hi 


PPjVjTl 1 ) < 


H{Vp) 


\D 2 h\\L°°\v \ 2 < 


C\\P\\] 


This provides the desired estimate on the boundary term. 


□ 


Remark A.5. In the above proof, we needed to use the L°° norm of p only in the boundary term. 
When there is no boundary term, the L 2 norm is enough, in order to handle the term (// a ). In 
both cases, the norm of p can be bounded in terms of the initial norm multiplied by e ct , where 
C bounds the divergence of v. On the other hand, in the torus case, one only needs to suppose 
Po £ L 2 and in the convex case po £ L°°. Both assumptions are satisfied in the applications to 
crowd motion with density constraints. 

We have seen that the constants in the above inequality depend on e and explode as e —>• 0. This 
prevents us to obtain a clean estimate on the BV norm in this context, but at least proves that 
Po £ BV pt £ BV for all t > 0 (to achieve this result, we just need to take e = 1). Unfortunately, 

the quantity which is estimated is not the BV norm, but the integral / iL(V p). This is not enough 

Jn 

for the purpose of the applications to Section 4, as it is unfortunately not true that the projection 
operator decreases the value of this other functional 2 . 

If we stay interested to the value of the BV norm, we can provide the following estimate. 

Lemma A.6 . Under the assumptions of Lemma A. 4, if we suppose po € BV(Ll) n L°°(U), then, 
for t <T, we have 

(A.4) f |Vpf|dx< f |Vp 0 | dx + CVt, 

J Q J Q 

where the constant C depends on v, on T and on H/OqIIl 00 - 


2 Here is a simple counter-example: consider fi = g{x) da: a BV density on [0, 2] C R, with g defined as follows. 
Divide the interval [0,2] into 2 K intervals Ji of length 2r (with 2 rK = 1); call ti the center of each interval A (i.e. 
ti = i2r + r, for i = 0,..., 2K — 1) and set g(x) = L + \Jr 2 — (x — ti) 2 on each J; with i odd, and g(x) = 0 on Ji for 

i even, taking L = 1 — 7rr/4. It is not difficult to check that the projection of fi is equal to the indicator function of 

the union of all the intervals Ji with i odd, and that the value of f H(Vp) has increased by K(2 — n/2)r = 1 — tt/4, 
i.e. by a positive constant (see Figure 3). 
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Figure 3. The counter-example to the decay of 
the total legth of the graph 



, which corresponds to 


Proof. Using the L°° estimate of Lemma 4.2, we will assume that H/o*|jis bounded by a constant 
(which depends on v, on T and on HpoIIl 00 )- Then, we can write 

[\Vp t \dx< ( H(\7p t )dx< f H(Vpo) dx + Cet + — < [ (\S7p 0 \ + e) dx + Cet +—. 

Jn J n Jn £ Jn £ 

It is sufficient to choose, for fixed t, e = y/t, in order to prove the claim. □ 

Unfortunately, this y/t behavior is not suitable to be iterated, and the above estimate is useless 
for the sake of Section 4. The existence of an estimate (for v Lipschitz) of the form TV(pt) < 
TV(po) + Ct , or TV(pt) < TV(po)e Ct , or even f(TV(pt )) < f(TV(po))e Ct , for any increasing 
function / : R + —> M + , seems to be an open question. 
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